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Solving the Infinite-horizon Constrained LQR 
Problem using Accelerated Dual Proximal 

Methods 

Giorgos Stathopoulos, Milan Korda, and Colin N. Jones. 


Abstract 

This work presents an algorithmic scheme for solving the infinite-time constrained linear quadratic 
regulation problem. We employ an accelerated version of a popular proximal gradient scheme, commonly 
known as the Forward-Backward Splitting (FBS), and prove its convergence to the optimal solution 
in our infinite-dimensional setting. Each iteration of the algorithm requires only finite memory, is 
computationally cheap, and makes no use of terminal invariant sets; hence, the algorithm can be applied 
to systems of very large dimensions. The acceleration brings in ‘optimal’ convergence rates 0(l/fc^) for 
function values and 0{l/k) for primal iterates and renders the proposed method a practical alternative 
to model predictive control schemes for setpoint tracking. In addition, for the case when the true system 
is subject to disturbances or modelling errors, we propose an efficient warm-starting procedure, which 
significantly reduces the number of iterations when the algorithm is applied in closed-loop. Numerical 
examples demonstrate the approach. 


Index Terms 

Constrained LQR, Alternating minimization. Operator splitting 


1. Introduction 

An important extension of the famous result of [[T| on the elosed form solution of the infinite- 
horizon linear quadratie regulation (LQR) problem is the case where input and state variables 
are constrained. This problem is computationally significantly more difficult and has been by 
and large addressed only approximately. A prime example of an approximation scheme is model 
predictive control (MPC), which approximates the infinite-time constrained problem by a finite¬ 
time one. Stability of such MPC controllers is then typically enforced by adding a suitable 
terminal constraint and a terminal penalty. The inclusion of a terminal constraint limits the 
feasible region of the MPC controller, and, consequently, the region of attraction of the closed- 
loop system. In practical applications, this problem is typically overcome by simply choosing a 
“sufficiently” long horizon based on process insight (e.g., dominant time constant). Closed-loop 
behavior is then analyzed a posteriori, for instance by exhaustive simulation. 

There have been few results addressing directly the infinite-horizon constrained LQR (CLQR) 
problem. Among the most well-known efforts are the works Q, Q and Q. The authors of 
Q suggest a scheme for offline computation of a sufficient horizon length. The solution of 
the corresponding quadratic program (QP) is then equivalent to the original infinite dimensional 
problem. The reported results are somewhat conservative, while the offline part of the proposed 
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algorithm can be computationally prohibitive since it involves the solution of a possibly non- 
convex problem, the computation of a positively invariant set as well as a vertex enumeration 
problem. The authors of Q extend the work of Q by solving a sequence of QPs of finite 
horizon length, which is monotonically non-decreasing. After each QP has been solved, the 
inclusion of the final state in a positively invariant set associated to the optimal unconstrained 
LQ controller is checked; if the final state is not included in the set, the horizon was insufficient 
and has to be increased. Finally, the authors in Q employ parametric quadratic programming and 
a reachability analysis approach to compute the least conservative horizon length that ensures 
optimal infinite horizon performance. Although this work provides the exact necessary horizon 
length for the feasible set of initial states, it suffers from tractability issues since it is based on 
state-space partitioning and thus can only be applied to small systems. 

Our approach is inspired by the framework of proximal gradient methods, a class of algorithms 
known for their simplicity when it comes to solving convex optimization problems with simple 
constraints [|^. From this family of algorithms, we use the (accelerated) version of the Forward- 
Backward Splitting method ((A)FBS) (see, e.g. [|7] Chapter 25] and Q). The idea is to condense 
the problem and describe it in terms of its input variables, then formulate the dual of the corre¬ 
sponding (infinite-dimensional) QP, and apply the (A)FBS method to solve it. More specifically, 
the method decomposes the QP into two subproblems, the first one being an infinite-dimensional 
least squares problem and the second one a simple clipping of an infinite sequence to the non¬ 
positive orthant. The subproblems are solved repeatedly (with the solution of one influencing 
the cost function of the other) until convergence to the solution of the original problem. This 
is in contrast to the approach of [[^, which requires the solution of a sequence of constrained 
QPs. We show that both sub-problems of the proposed algorithm can be solved tractably (which 
is not a priori obvious since we are working with infinite sequences), the first one by solving a 
finite-dimensional system of linear equations (with the possibility to pre-factorize the matrices) 
and the second one by simple clipping of finitely many real numbers on the non-positive real 
line. Convergence of the scheme (with rate 0{l/k‘^) for function values and 0(l/k) for primal 
iterates) to the optimal infinite-horizon sequence is guaranteed under mild assumptions. Therefore 
the proposed algorithmic scheme provides a means to compute the solution of the infinite-horizon 
constrained LQR problem with guaranteed convergence. 

This work is based on our recent result in Q, where the same problem is solved by employing 
a primal-dual scheme, the Alternating Minimization Algorithm (AMA) [10|, which is equivalent 
to FBS. The approach is equivalent to splitting the infinite-horizon constrained LQR problem into 
an unconstrained LQR problem and a proximal minimization problem. In the current work, the 
result is developed in both theoretical and practical directions. From the theoretical perspective, 
we (i) drop the limiting assumption for positive-definiteness of the state penalty matrix Q that 
existed in [[^, and, more importantly, we (ii) prove convergence of an accelerated counterpart of 
the original FBS method, a fact that enables a (significantly) faster convergence rate. From the 
practical perspective, we provide a fully implementable method, competitive for real-time control. 
We (i) eliminate the need for knowledge of uncomputable quantities, (ii) propose computationally 
efficient ways to solve the optimization subproblems and (hi) propose a warm-starting scheme 
that performs well in practice. 

The paper is organized as follows: In Section |I^ we introduce the problem and express it 
in terms of its dual variables by means of the proximal splitting framework. In Section 1111] we 


present the accelerated forward backward splitting algorithm to solve the problem and show 
that each iteration of the algorithm can be carried out tractably. The convergence proofs for this 
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scheme are given in Seetion while Seetion |V] diseusses the eomputational aspeets that make 
the algorithm praetieal to use. Seetion VI presents two numerieal examples: A toy example 
of an unstable system with two states and one eontrol input illustrates the main features of 
the algorithm. Subsequently, we demonstrate the praetieal applieability of the algorithm on a 
linearized model of a quadeopter with 12 states, 4 inputs and polytopie eonstraints. Finally, 
Appendiees and provide the proofs for the results presented in Seetion |W 


II. Problem statement and dualization 
The eonstrained regulation problem ean be written in the general form 


minimize 


i=0 


subjeet to Xj+i = Axi + Bui, 

To Tinit 
CxXi ^ Cx 
CuUi ^ Cu , 


i e N 


( 1 ) 


with variables Xi G and Ui G M"*, and data Cx G and c„ G 

We make the following standing assumptions: 

Assumption 1: The pair {A, B) is stabilizable, the optimal value of problem Q is finite, the 
set 

A := {x G M"" I CxX < Cx} 

eontains the origin in its interior, the matrix has full eolumn rank, the matrix Q is positive 
semidefinite and R is positive definite. 

Remark 1 (Stability): Clearly, under Assumption if the problem Q is feasible and the pair 
(A, y/Q) deteetable, then the eontrol sequenee optimal in Q is stabilizing. Therefore, there is 
no need to enforee stability in an ad hoe way as is eommonly done when the infinite-time 
problem Q is approximated by a finite-time one solved in a reeeding horizon fashion. 

From now on we write infinite sequenees and infinite-dimensional operators in bold font. 

The problem ean be rewritten in the dense form, i.e., by writing the states as funetions of the 
inputs. This is done by defining the operators 


A = 


Q = diag(Q,Q,...) , R = disLg{R, R,...) , 

Cx diag(C*j;, Cxy • • •) ) Cx)...) > 

C,= 


r A ~\ 


■ B 

0 

0 ■ 


A 

A 2 


AB 

B 

0 ■ 



, B = 

A^B 

AB 

B ■ 


L J 










. [C.S]. . 




H = B*QB + R, G = A*QB, F = A*QA + Q , 
C=[C* C* ...]*, c=[ct,c*,...r , 


( 2 ) 


where we denote by e* the (infinite dimensional) row veetor with only one nonzero element 
equal to one at position i, by [■] • the bloek row of size pa: x cx) of the eorresponding operator 
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and by * the adjoint of an operator (i.e., the infinite-dimensional analogue of transpose; see 
Appendix A for a brief introduction to operators). Using the above, 0 can be written in the 
form 

minimize ^u*Hu + h*u + r ,, 

subject to Cu < c , 

where h = G*Xiniu r = u is the infinite sequence u := [ mq , .. .]* and H: T-Lu 

V-u, C: T-Lx —)■ T-Lu, c G T-Lx, where Ws are suitable Hilbert spaces specified next: Any sequence 

z := • • •) 

is viewed as an element of an /^-weighted (or real Hilbert space (see Appendix A, 
Definition induced by the inner product 

OO 

{z,y) = '^w^zjyi , yye'H^,zen^, ( 4 ) 

i=0 

where w is an appropriately chosen weight (see Appendix A, Definition [^. The norm of any 
z G T-Lz is thus given by 


\z\\n. ■= = 


\ 


OO 


i=0 


W'-\\Zi 


Unless stated otherwise, for the rest of the paper by a Hilbert space we mean the real Hilbert 
space as just introduced. 

The dual of problem Q can be derived by making use of Legendre-Fenchel duality (see, e.g.. 
Chapter 7, GD)' By defining 

f{u) = -u*Hu + h*u + r, g{Cu — c) = 6-{Cu — c) (5) 

and (5_(-) being the indicator function for the nonpositive orthant, we can rewrite Q as 

min {f{u) + g{Cu - c)} 


and then express the Lagrange dual problem by using the dual variable A as 


max < min {f{u) + (c — Cu, A)} — 


A L u 




max 

A 


max{(C*A,u) - f{u)} - g{X) + (A,c) \ ^ 


max{-/*(C*A)-^(A) + (A,c)} , 

The solution can be derived by solving 

imn{r{C*\)+g{X)-{X,c)} . 

A 

The function involved in the minimization problem Q is 

/*(C*A) = ^ (C*A, H-^C*X) - (C*A, H-^h) + ^ {h, H-^h) - r 


( 6 ) 

(7) 
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Finally, the problem we are interested in solving ean be east in a more eompaet form as: 

minimize F{\) := h*{\) + (8) 

with variables A, h*{\) = /*(C*A) — (A, c) and /i*(A), (5_(A) being proper lower semi- 
eontinuous eonvex funetions from T-ix to M U {+oo}. 

Before proeeeding, we elaborate on the reasons why the original problem Q had to be 
reformulated in order to solve it. There are two reformulations, namely, posing the problem as a 
funetion of the input sequenees only, resulting in Q, and dualization of Q, resulting in ([^. The 
reason for eonsidering the eondensed formulation is the need for strong convexity of the primal 
objeetive, whieh implies the Lipsehitz eontinuity of 0 Corollary 18.16]. By using the 

eondensed form we avoid the restrietive assumption of Q 0 required in Q. 

The reason for eonsidering the dual problem is simplieity in the evaluation of the proximal 
operator (see Appendix [a} Defintion|^ of the funetion ^-(■), whieh is a simple projeetion on 
the non-positive orthant (i.e., eomponentwise elipping) as opposed to the primal ease where one 
would have to projeet on a generie polytope of the form Cu < c. 

Both steps are erueial for sueeessfully applying the Forward-Backward Splitting (FBS) method, 
as presented and analyzed in 0 Theorem 25.8] 


III. Solution using AFBS 

Problem Q is a eomposite minimization problem (i.e., a minimization of a sum of a smooth 
and a non-smooth funetion) and will be a solved using an accelerated forward-backward splitting 
(AFBS) method, whieh is an aeeelerated variant of the forward-baekward splitting proposed in 0 
Theorem 25.8]. The aeeeleration eomes from a Nesterov-like momentum sequenee from p^ . 

This modifieation allows aeeeleration of the FBS method in a very simple manner and adds 
praetieally zero eomputational eomplexity per iteration, in eomparison to the original version. 


The idea ean be traeed baek to Polyak and the so ealled heavy ball methods [ 13| for minimizing 
a smooth eonvex funetion /(■): 

= A^ + a^(A^ - 
Afc+i _ A^ _pfcv/(A^) , 

where G [0,1) is an extrapolation faetor and a stepsize parameter. This seemingly small 
ehange of updating the new iterate as a linear eombination of the two previous iterates greatly 
improves the performanee of the original gradient seheme. 


In his seminal paper [14|, Nesterov modified the heavy ball method by simply evaluating the 
gradient at the extrapolated point A^ instead of A^. In addition, he proposed a speeial formula 
for eomputing the relaxation sequenee resulting in an optimal eonvergenee rate for 

the seheme. Subsequently, Giiler extended Nesterov’s results to the proximal point algorithm. 


handling the minimization of the sum of a smooth and a nonsmooth funetion [15|. 

The aeeeleration that optimal first order methods enjoy eomes solely from the additional 
momentum of the minimizer sequenee due to the relaxation sequenee Among the many works 
that derived optimal relaxation sequenees similar to that of Nesterov, a distinguished one is that 
of Beek and Teboulle in who used Nesterov’s framework to aeeelerate the FBS method. 
The resulting algorithm, eommonly known as FISTA (Fast Iterative Shrinkage Thresholding 
Algorithm) gives an optimal, eonvergenee rate in terms of the funetion values. This result 
has valuable praetieal and theoretieal implieations, however, eonvergenee of the iterates of the 
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FISTA algorithm has not yet been proven. Building upon FISTA, Lorenz and Poek suggest in [[^ 
a modified version of the method that aehieves weak eonvergenee of the iterates to the solution of 
the optimization problem in a Hilbert spaee setting, yet not proving anything about eonvergenee 
of the funetion values. 

The several ‘optimal’ relaxation sequenees have been put under a eommon framework in 
the reeent work [12|. The authors showed that any sequenee of the form 

with a > 2 satisfies the inequality k > 2. Then the sequenee defined as 

allows for the optimal 0{l/k‘^) eonvergenee rate in terms of the funetion values, as 


a = 


tk+i 


well as weak eonvergenee of the iterates, and henee provides an optimal first order seheme. We 
denote the algorithm emanating from this seheme the Aeeelerated Forward-Baekward Splitting 
(AFBS) algorithm, and write it down for problem @ as follows: 


Algorithm 1 AFBS for Problem (8) 


0: Initialize A° = 0, a > 2, = 0, 

L(h* *) — a Lipschitz constant of V/i*. 

repeat 

2: + a^(A'= - A^"^) 

3: A"+' = min{A' - 0} 

until termination eondition is satisfied. 


The iterative seheme above is very simple: in this ease the AFBS boils down to a variant of 
the fast projeeted gradient method as proposed by Nesterov 02). In order to apply algorithm ([T]), 
we need to to be able to 

• evaluate the gradient of h*{-) 

• represent A^ and A using a finite amount of memory. 

The remaining steps of the algorithm are simple sealar or veetor updates or eomponentwise 
elipping on the non-positive orthant (Step 3), both of whieh ean be earned out inexpensively 
provided that A^ ean be represented using finite amount of memory. 

In the rest of the text -Plq and A^lq will denote the positive-semidefinite solution to the 
diserete-time algebraie Rieeati equation and the eorresponding LQ optimal state feedbaek matrix 
assoeiated with the matriees {A, B,Q, R). 

We start by finding the gradient of To this end, define the Lagrangian of the (infinite 
dimensional) problem Q, written in its equivalent dense form and truneated at any T > 0, as 


C{u, X, A|T) 


T-l 

XtPlqXt + '^^J Q^i + uj Rui - te* 

i=0 


CxXi Cx 
CqjUi C71 


T 


A*. 


Lemma 1: If A is sueh that A* = 0 for all i > T, then 

CxXi Ct 


[Vh*{X)], = 




where 


{x, ii) = argmin{£(u, x, A|T) | Xj+i = Axi + Bui, Xq = Xinit} 

U^X 


(9) 

( 10 ) 
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for i e {0,..., T}, and 


Xi = {A + BKi^qY '^Xt, Ui = Ki^qXi 

for i > T. 

Proof: This is a standard result from duality (see, e.g., Chapter 5]), notieing that 
the minimization in ( [TO] ) is equivalent to the minimization of the Lagrangian of Q under the 
assumption that A* = 0 for all i >T. ■ 

Lemma gives us a way to eompute the gradient of h*. Clearly, this gradient is an infinite 
sequenee and therefore eannot be stored direetly, but it is available to us explieitly for i G 
{0,... T} and implieitly, through the dynamies of the system x+ = {A + BKi^q)x, for i > T. 

Now we show that and A ean be represented using a finite amount of memory. 

Lemma 2: If A^ and A are generated by Algorithm then for eaeh k there exists a < cx) 
sueh that Af = 0 and A^ = 0 for all i >T^. 

Proof: We have A° = 0 and A° = 0 for all 7 > 0 and henee = 0. Assume now fc > 0 
and A^ = 0 and Af = 0 for all i > T^. Then, evaluating V/i*(A^) using Lemma with T = T^, 
we see that the sequenee {xi,Ui) is generated by the uneonstrained LQ eontroller for i > 
and henee eonverges to the origin. Sinee the set X has the origin in the interior we eonelude 
that there exists a time < oo sueh that C^Xi < Cx and CuUi < Cu for all i > We 

define = max{T^, In view of (|9), we eonelude that A — £^^V/i*(A ) > 0 for 

all i > Therefore, in view of Step 3 of Algorithm we have = 0 for all i > 

Finally, A ^ is a linear eombination of A^ and A^^^ and henee = 0 for all i > T^+^. ■ 
To determine eomputationally (given and and u~^^) we simply find the first 


time T'^ that 


enters a given subset S, with 0 G int S, of the maximum positive 


CuK 

Cx 


X < 


Cu 

C-x 


y invariant 
. The time 


u. 


k+l 


< Cu 


set of the system x+ = {A + BKi^q)x subjeet to the eonstraint 

is then equal to the first time no less than sueh that CxX^llli < Cx and C, 
simultaneously hold for all i G ... ,T^}. More formally, we have the equality 

= min {T>T’^\3T^ s.t. CxX^+^ < c^, Cuu’l+^ < c„ Vi G {T, ..., T^}, and 4+^ G S}. 

( 11 ) 


Remark 2 (Computation ofT^): In praetiee, to determine after solving ( |T^ , we iterate 
forward the system dynamies x+ = (A + BKi^q)x starting from the initial eondition until 

4+^ G 5. 

Remark 3 (Set S): The set S is determined offline and is not required to be invariant. A gooc 

eandidate is the set {x | x~^P i,qx < 1} sealed sueh that it is ineluded in |x 

or any subset of this set eontaining the origin in the interior. 

Now we are ready to formulate an implementable version of the abstraet algorithm 


CuK 

a 


X < 


C-n 

Cfr 
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Algorithm 2 AFBS for the CLQR problem 


Require: xinit, Q ^ 0, i? >- 0, of full eolumn rank, w = min |l, ^2 ^ | 


, a > 2, Q!° = 0 


a: Determine -Plq, A'lq solving the uneonstrained LQR 
problem assoeiated with the matriees (A, B, Q, R). 
b: Determine a set S, with 0 G int S, ineluded in 
any positively invariant set for the system 
= (A + BKi^q)x subjeet to the eonstraint 


CuK 

a 


X < 


Cqj 

Cr 


. See Remark 


e: Initialize A'^ = 0, = 0, a > 2, 

> 0 or L{h*) — a Lipschitz constant ofVh* (optional). 
for fc = 0,... do 
l:a'= = |^, k>l 
2: Af = Af + a^(Af - i = 1, 

3: Set 


Rk 


[x 


k +1 - fe+1 


u ' -) = argmin{£(u, x, A \T^) \ Xj+i = Ax* + Bui, xq = Xinit}, i = 0,..., T'* 

U,X 

i*p = {A + BirLQ)4‘+‘. »> r‘ 


4: Determine (see Remark 
6: Choose stepsize p (see Remark 

7: Set A,^+^ = min ( A,^ - 


w 




CuUj^ 


k +1 


,0 


8: If termination eondition is satisfied, solve KKT system (see Remark 

end for 


Remark 4 (Stepsize): In Step 6 of the algorithm, a stepsize is seleeted. One option is to fix 


a eonstant stepsize p = 


L{h*y 


whieh needs a global Lipsehitz eonstant of V/i*(-). This ean 


be eomputed offline (see Seetion V-A). Alternatively, one ean use a backtracking stepsize rule, 


inspired from p^ , whieh ean be used in eombination with the global estimate. The proeedure 
is analyzed in Seetion |V| 

Remark 5 (Role of the weight w): It is worth mentioning that working in the weighted Hilbert 
spaee is mueh more than a mathematieal formalism and has serious praetieal implieations. 
In the ease of unstable systems, a nontrivial sequenee of weights has to be ehosen sueh that 
the growth of the largest unstable eigenvalue of the state matrix A is bounded by a faster 
deeaying sequenee so that the operator C remains bounded, an assumption neeessary for applying 
the proposed method. At the same time the sequenee of weights will aet as a left and right 
preeonditioner on the Hessian operator of the quadratie form in ([^, as well as a left preeonditioner 
of C. This sealing ean seriously affeet the numerieal performanee of the proposed algorithms, 
as we will see in subsequent seetions. Note that for stable systems the sequenee ean be trivially 
set to 1, henee no sealing oeeurs. These elaims are explained in more detail in Appendix]^ 

Remark 6 (Termination): Algorithm terminates when a prespeeified aeeuraey is reaehed in 
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terms of the progress of the dual sequence. The extracted primal sequence is given by 

(a;^, u^) = argmin{£(w, x, A^|T*^) | Xi+i = Axi + Bui, Xq = Xinit} (12) 

U,X 

for i G {0,T^} and x^ = {A + x^k for i > T^. In Theorem below it is proven 

that (a;^, u^) tends to the optimal constrained LQ solution. At any finite iterate, however, the 
sequence (a;^, u^) may violate the constraints. In order to remedy this we solve upon termination 
an equality-constrained QP where we minimize the objective function subject to the active 
constraints at optimality. The active constraints can be (approximately) detected by looking 
at the nonzero values of the dual vector at termination. This step comes at a very small cost 
since it involves one solution of a KKT system of linear equations. 


IV. Convergence results 

In the previous section we gave an implementable algorithmic scheme that computes the 
solution to the CLQR problem. Here we provide all the necessary proofs which allow us to 
assert that the solution generated by Algorithm via ( [T^ indeed converges to the true optimizer 
of the CLQR problem. In what follows A°° denotes any optimal solution to the dual problem ([^ 
(which exists under Assumption but may not be unique) and the optimal solution 

to the primal problem ([T]). Our main result is: 

Theorem 1 (Main Theorem): Suppose Assumption holds and let A^ be a sequence of iterates 
generated by Algorithmand (x^^u^) the associated primal sequence given by (12) and let L 
be a Lipschitz constant of V/i*(-). The following statements hold: 

(i) The composite function F{\) = h*{X) + ^-(A) as defined in ([^ converges as 

(ii) The sequence of the dual iterates converges weakly (see Definition in Ap¬ 

pendix 1^ to an optimizer, that is, 

A^ ^ A“ 

for some A°° G argmin(F). 

(iii) The input sequence converges strongly to the unique minimizer as 




u 


m-u. ^ 


< a 



mx 


(k + a — 1) 


where /x > 0 is the strong convexity modulus of f{u). 

(iv) The state sequence converges strongly to the unique minimizer as 


\x’^-x°°\ 




< a\ 


|B||2L||A°-A" 


Wx 


/X (Jx, + a — 1) 


(v) The sequence is bounded. 

Proof: 

(i) Convergence of F{X^) with a constant stepsize is proven in [12 Theorem 1]. Convergence 
at the same rate with an adaptive stepsize generated from the backtracking Algorithm is 
proven in Lemma Appendix [D| 
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(ii) The proof is stated in [12 Theorem 3]. 

(iii) The idea is to upper bound the input sequenee’s eonvergenee rate making use of the result 
in point (ii). In order to do so we make use of strong duality. The proof is inspired from 
[ [T^ Theorem 4.1] and is as follows: 

Let < 0 generated from Step 7 of Algorithm Denote 

= argmin {/'(u) :=/(li) + (A^, c — Cu)} , (13) 

where f{u) = }^u*Hu + h*u + r as defined in Seetion|^ Then we have that the Lagrangian 
of ([^ evaluated at A^ is C{u,X^) = f'{u). The funetion f{u) is strongly eonvex with 
modulus /i > Amin(-R) > 0, where Amin(-R) denotes the smallest eigenvalue of R. Strong 
eonvexity of f'{u) with modulus /i follows direetly. Using ([ Ts ]), it holds that 




u 


fc ||2 




Wu G Rr 


or, equivalently, 


Ciu,X’^) - > ^\\u 


u 


k \\2 


m-u.-' 


'iu G 


(14) 


Substituting u = u°° in (14) and by observing that max C{u°°,X) > C{u°°,X^), we have 
that 


Ciu^,X^)- C{u\X^) > ^\\u 


u 


k \\2 




Vw G n, 


(15) 


We have managed to derive an upper bound for the distanee of the generated sequenee 
of primal minimizers from the optimal one. The last step is to show that the 

Lagrangian C{u,X) is assoeiated to the eomposite objeetive F{X). This ean be easily 
shown as follows: 


A^) = min {f{u) + (A^, c — Cu) } 

= — max [—f{u) + (A^, Cu) \ + (A^, c) 

= -/*(C*A^) + (A^c) 

= -F(A'^), by (§ . 


From strong duality and the faet that —F(A^) eonverges to the optimal dual value (point 
(i)), we have that the optimal value of the dual funetion —F{X) eoineides with that of the La¬ 
grangian evaluated at the saddle point {u°°, A°°), i.e., C{u°°, X°°) = max {—F(A)} = —F{X° 

_ _ Xg'Ha 

(see Seetion 5.5.5]). Making use of point (i), inequality ([15|) beeomes 




11 ftz, 


< F(A^) - F(A“) < 


a^LllA" 


\CXD ||2 

WHx 


2{k + a-iy 


(16) 


whieh eoneludes the proof. 

(iv) The state sequenee is generated by 


= Ax\ 


init 


Bu'^ 


(17) 
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Strong convergence of the input sequence along with the facts that B: Hx —^ 'Hu 

is bounded (follows directly from Lemma in Appendix C) and the uniqueness of u°° 
prove strong convergence of the state sequence with rate 1/k, i.e., 

< \\B\\\\u'^ - 


< a\\B 



Wx 


/r (fc + a — 1) 


with the last inequality following directly from point (iii). 

(v) The proof is already presented in Q, but we repeat it for completeness. First note that for 
the statement to hold it is sufficient to show that 


limsupT'' < oo. 

k^oo 


(18) 


To prove (18), define the sequence of the first hitting times of the interior of S as 

:= inf{i > | E int 5}, keNU {+cx)}, 

where r°° < cxd is the hitting time of the optimal state sequence x°°. Clearly, 
and H < oo since the origin is in the interior of S and for each k E N the sequence 
(a;f)igN generated by Algorithm converges to the origin as i ^ oo. We shall prove that 
limsup^_^o^ <T°° < oo, which implies (18). 

The key piece for the proof is the weak convergence of x^ to x°°. This claim is proven in 
Lemma Appendix and is a consequence of other intermediate results, all presented in 
the same Appendix. For the purpose of contradiction assume that there exists a subsequence 
i £ with limj_).oo + 1. Since the sequence of hitting times is integer 

valued, this implies that there exists a j* G N such that + 1 for all j > j*. We 

now use this to contradict the weak convergence of x^ to x°°. To this end, observe that 
X E int S whereas f ^ int S for all J > f- By the definition of the interior there exists 


an e > 0 such that z E int S for all 2 ; with 
for all j > j*, and consequently 


\z — x^ 


I 2 < e. Therefore I lx.,. 


— x: 


2 ^ 


> e 


— X°°,z) = w^(x^^ 


i=0 


> {x% — X^oo ) ' 

> > 0 , 
k 

for a sequence z with Zr<=° = X ./00 — x!^ and for all j > j*, contradicting the weak 
convergence of x^ to x°° asserted by Lemma 


V. Computational aspects and warm-starting 

Having presented the algorithm and its convergence results, we now focus on the computational 
aspects that render the algorithm a practical option to alternatives such as MFC. We start with 
explaining how no prior knowledge of a fixed stepsize is needed, give some references concerning 
the solution of the linear system, which is the most expensive operation of the method, and we 
conclude the section by suggesting a warm-starting scheme. 
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A. Stepsize selection 

The stepsize used in Algorithm is eomputed as the reciproeal of the Lipsehitz constant of 
h*. For the problem discussed here this can be explicitly computed. We have that 

||V/i*(Ai) - V/i*(A2)|| = \\CH-^C*W{\i - A2)|| 

< ||Cif-^C*iy||||Ai - A2II 

= ||if-i||C|n|||Ai-A2|| , 


where W is the diagonal operator constituted of the decaying weighting sequence, i.e., W = 
diag(/, teJ, ...). The last equality follows from the fact that W contains a non-increasing 
sequence with the largest element being one. Hence L{h*) = ||iT~^|| HC'lp, which requires com¬ 
putation of the operator norm ||iT~^|| and C. The proofs for boundedness of the two operators, 
as well as the computations of their bounds are derived in Appendix Lemmas and |7} 
Although valid, this offline computation of the stepsize tends to be conservative in many cases, 
due to the conservativeness of the computed upper bounds. An elegant and practical method to 
achieve faster convergence is to employ an algorithm that locally estimates the Lipsehitz constant 
online, at every iteration of Algorithm]^ In order to do so, we use the backtracking stepsize rule 
suggested in p^ . The idea is simple: after each iteration of the algorithm, we make a quadratic 
approximation model of the function around the successor point, making use of the knowledge 
of the exact point and its gradient value. A quadratic term with varying curvature is added on 
top of the linear (first order Taylor) approximation, and the curvature is adapted recursively until 
our quadratic approximant upper bounds the original function, centered around the given point. 
Thus, the quadratic model’s curvature is an estimate of the Lipsehitz constant of the gradient 
of the original function. It is proven in [16| that the locally evaluated Lipsehitz constant is 
related to the global one by l3L{h*) < where (3 = and 7 > 1, > 0 being 


an initial estimate. Consequently the rule allows for smaller L’s and hence larger stepsizes, i.e., 
faster convergence in practice. 


B. Complexity 

The most expensive operation of Algorithm is the linear system solve in Step 3. There is a 
variety of ways to perform this step, i.e., solving the KKT system or perform the Riccati recursion 
when both states and inputs are considered, invert the dense Hessian when only the inputs are 
considered. In the first case, a sparse (permuted) LDL^ factorization can be performed with cost 
T{n-ymY flops, followed by a forward-backward solve at T(n -f m)^ flops. This approach has 
been followed i n p0| . A discussion on the KKT system solve and the Riccati recursion approach 
is contained in pi[|, where the corresponding complexities are analyzed and compared in detail. 

In the case of the condensed formulation, the linear system solve can be efficiently perfomed by 
first applying a Cholesky factorization on H (being a finite truncation of the H operator in Q), 
followed by a forward-backward substitution, (see Appendix C]). Although the condensed 
form of the optimal control problem that is used in the derivations is, generally, not advised 
for long horizons, recent advancements can render this approach very efficient [22|, [231. More 
specifically, the two proposed algorithms that perform factorization and solve of the condensed 
system in [231 come with a reduced complexity of 0{Trc‘) and O(T^n^), respectively. 

Whether considering the sparse or the dense formulation, note that the factorization steps 
would have to be performed several times until the ‘correct’ horizon has been identified, 
since the size of the corresponding matrices (KKT or Hessian H) increase as the algorithm 
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progresses. This typieally happens within the first few tens of iterations. An alternative would 
be to direetly perform a matrix inversion (for small to medium sized systems) and only faetorize 
onee, when the horizon seems to have stabilised. 

The baektraeking seheme eontributes to the eomplexity of Algorithm by requiring a number 
of funetion evaluations per iteration, both for the quadratie model and for the original smooth 
funetion h*{-) (see [16| for more details). The rest of the steps are simple veetor updates of 
negligible eost. 


C. Warm-starting 

In the nominal ease, i.e., when no noise and no model uneertainty are present, the open 
loop infinite-horizon eontrol sequenee generated from Algorithm eoineides with the eontrol 
sequenee generated by the optimal elosed-loop feedbaek eontroller. Consequently, there is no 
need to re-optimize in a reeeding horizon fashion. Solving the CLQR problem just onee is 
suffieient. 

In the more realistie seenario where the predieted initial state differs from the measured one, 
the algorithm has to be re-applied. Provided that the predietion is not very different from the 
aetual state, a good strategy is to initialize the deeision variables (states and inputs) of the 
new problem with the values predieted from the previous one. This is eommonly known as 
warm-starting. In our ease, warm-starting has to be performed in the dual variables. 

Note that onee Algorithm has run onee, a hitting time T°° has been generated, along with 
an optimal dual sequenee of eorresponding length. It is expeeted that when eomputing the 
eontrol law the hitting time should deerease by one at eaeh solve, provided warm-starting from 
the optimal dual sequenee that was generated onee in the beginning. Henee, we suggest a heuristie 
seheme where the ‘eonstrained’ (nonzero) part of the preeeding shifted dual sequenee is used 
to initialize eaeh subsequent CLQR problem. The eomputation time thus reduees signifieantly, 
with the horizon praetieally shrinking to zero onee the intial state is identified to be inside the 
maximal positively invariant set of the LQ eontroller. 

Applieation of the seheme is presented in Seetion |V^ It is observed that it behaves partieularly 
well for small perturbations of the initial state. 

VI. Examples 

For illustrative purposes, we run the algorithm on two systems, a small unstable system with 
two states and one input and a linearized model of a quadroeopter with 12 states and 4 inputs. We 
use the small example as a benehmark for graphieal illustrations, while the larger one exhibits the 
eomputational effieieney of the proposed seheme. The eomparison is performed against the same 
implementation of the AFBS algorithm for finite horizon lengths. It is of eourse understood that 
there exist several methods eapable of solving a finite horizon MFC problems (see interior point, 
aetive set, ete.), among whieh, optimal first order methods have gained eonsiderable attention over 
the last few years, rendering them a eompetitive alternative p4| , p5| . Consequently, eomparing 
against an optimal first order method provides a valid basis for evaluating the potential of our 

seheme. In the two examples the termination eriterion is simply set as ||A^ — A || < 10“^. 

A. Toy system 

Consider the following system defined as 

O' 

0.0787 ’ 


A = 


1.1 

0 


2 

0.95 


B = 
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Hitting times T" 

Fig. 1. Histogram of T°° = maxfc{T^} for 750 initial conditions of the 2 state system sampled from a normal 
distribution around (—3,0.3) with covariance matrix diag(4, 0.4). 


Xi+i = Axi + Bui, 


with constraints 


^||oo — Halloo — 1 


and Q 



21. Note that the system is unstable, henee a nonzero sequenee 


of weights has to be ehosen in order to ensure boundedness of the C operator. We ehoose 
w = 1/1.1^ and the value of a in Step 1 of Algorithm is set to 4. 

The system is simulated for 750 different initial eonditions xq. In Figure the distribution of 
T°° = maxfc{T^} is depieted. We see that goes up to 30, while the mean value is 9. 

We eompare our method to the MFC approaeh both from the eontrol and the algorithmie 
performanee perspeetive. Regarding the former, we perform a eomparison of the feasible sets 
for finite horizon implementations versus the maximal eontrol invariant set in the ease of CLQR, 
as well as the optimal value of the eost funetion. Regarding the latter, we perform eomparisons 
in terms of the average number of iterations needed for eonvergenee for several horizon length 
values in MFC versus the CLQR ease. We also evaluate the eonservativeness of our approaeh by 
eomputing the aetual ‘optimal’ horizon length for eaeh initial eondition we simulate. For both 
CLQR and MFC we make use of the same AFBS algorithm with baektraeking employed and 
termination toleranee set to 10“^ as stated before. 

We perform the following simulations: we sample 243 initial eonditions Xi^, i = 1,... ,243 
from the maximal eontrol invariant set (31-step) of the aforementioned eonstrained system (see 
Figure [^. For eaeh point we eompute: 


1) The minimum horizon length sueh that x^^.^{xifi) resides in the maximum positively 
invariant set of the autonomous system x^ = {A + BKi^qjx, used as a terminal set. 

2) The hitting time T = T°° generated by our proposed seheme. 

The following seenarios are generated: Firstly, an MFC problem with terminal set and horizon 
length Tmin is solved, as deseribed above. Subsequently we remove the terminal eonstraint and 
solve the MFC problem again for the same horizon length T^in- We repeat the proeedure 
deseribed above (MFC with and without terminal set) for horizons T = 2Tmin and T = T*, 
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where T* is speeified by identifying the first horizon length for whieh the feasible solution of 
the MPC problem had the terminal eonstraint inaetive. Finally, a eomparison with the proposed 
CLQR approaeh is performed, with horizon T = T°°. In this way, six MPC problems with finite 
horizon as well as the CLQR problem are solved, for eaeh of the 243 initial eonditions. The 
total number of iterations is averaged by the number of eorresponding initial eonditions with the 
same minimum horizon length Tmin- The results are presented in Figure]^ 

The first observation from the plot is that inelusion of the terminal eonstraint generally 
inereases the number of required iterations sinee more eonstraints beeome aetive at optimality. 
This is espeeially the ease when the horizon is relatively short and the terminal eonstraint 
satisfaetion is imposed as we see in the blue eurves. As the horizon inereases (red and magenta 
eolor), the terminal eonstraint might be inaetive and there is no signiheant differenee in the 
number of iterations in the two oases. 

A trend of inorease in the number of iterations in all methods as T* inereases oan be observed. 
This is expeoted sinee inorease of T* amounts to sampling of the initial state from regions of 
the feasible set that are further away from the origin. A oonsequenoe of the latter is that more 
eonstraints beeome aetive at optimality. 

Figure 1^ depiets a generally expeeted behavior. Short horizons in oombination with a terminal 
set inorease signifioantly the iteration oount. As the horizon inereases the iteration eount drops. 
However, an unneoessarily long horizon (e.g., 2 Tru\u ) leads as well to an inorease. In that sense, 
one ean observe that MPC with T = T* behaves better than the other lengths. The eurve 
eorresponding to CLQR is denoted with blaok lines with asterisks. We observe that the method 
is eomparable to the MPC approaehes where T = Tmin and T = T* with terminal set, performing 
relatively better in small to medium horizon lengths and worse for larger horizon lengths. The 
reason for the latter is that the weighting operator W beoomes quiokly ill-oonditioned as T 
grows, sinee its diagonal elements deeay exponentially. This leads to ill-eonditioning of the finite¬ 
dimensional truneation of the dual problem ([^ sinee W preeonditions (through the weighted I 2 
inner produet Q) the objeetive of Q. As is eommonly known first-order methods struggle with 
ill-eonditioning. Effieient preeonditioning of the operators should improve this and is a topie of 
further investigation. 

We further perform the following simulation: The CLQR approaeh is applied, this time with 
the weights being set to one. The result is depieted as a blaek line with eireles. In this ease, 
CLQR elearly outperforms all MPC approaehes with terminal set, as well as most of those 
instanees without terminal set. Further simulations suggest that the proposed method performs 
very well, were we to drop the weights or when dealing with stable systems, when no sealing 
has to be performed. 

Another interesting point is that the average hitting times T°° generated from Algorithm 
are almost identieal to the averaged optimal ones T*, with a very slight inerease. This faet is 
depieted in Figure Q where the ratio is always very elose to one. On the eontrary, the minimum 
required horizon length is observed to be up to 45% smaller than the optimal length T* 
in some eases, leading to a signiheant inerease in the objeetive’s eost when eompared to the 
inhnite horizon approaeh. 

B. Quadcopter system 

The next system we eonsider is a quadeopter linearized in a hovering equilibrium. The system 
has 12 states whieh eorrespond to position, angle and the eorresponding veloeities. There are 
four inputs eorresponding to the four propellers. There are box eonstraints on all states and 


January 20, 2015 


DRAFT 


16 


Reach sets for given terminal set 



Fig. 2. Reachable sets for several horizon lengths and the LQ terminal set. The computation was done using 
the MPT3 toolbox | |2^ . It is apparent that a short horizon length reduces significantly the feasible region of the 
problem. 



Fig. 3. Comparison of MPC with finite horizon length, for several horizons, with CLQR. The x-axis corresponds to 
the optimal horizon length T* per initial condition. As T* increases, the corresponding states are sampled further 
from the origin. MPC with terminal set is depicted with the solid lines, while without terminal set with dashed 
lines. CLQR is performed both with the weight sequence, denoted as CLQR (provably convergent version), 
and without the weights, denoted as CLQR. 


inputs, mainly ensuring the validity of the linearized model. The system is marginally stable; 
thus the weight w is set to one. 

We simulate the algorithm starting from initial eonditions randomly seleeted as follows: 
starting from a random feasible initial condition, we generate random directions on a unit ball 
centered around it and sample points along each of them. The points are generated from a normal 
distribution with standard deviation 0.15. Finally, we keep the initial conditions that result in 
feasible closed Ioot problems. The result of this step is 272 feasible initial conditions for the 
CLQR Algorithm]^ A histogram of T°° = max^jT^} is presented in Figure 

We conclude the section by applying the warm-starting heuristic scheme suggested in Sec- 


Januaiy 20, 2015 


DRAFT 












17 


Horizons length comparison 



rri 'J~iOO _QQ 

Fig. 4. Evolution of the ratios -fr and for the sampled initial conditions. 


Hitting times distribution 



Hitting times T“ 

Fig. 5. Histogram of T°° = maxfejT^'} for 272 initial conditions sampled with a Hit-And-Run algorithm. 


tion |V| We consider two scenarios; in the first one we uniformly perturb the initial state by 
0.5% of its nominal value and run Algorithmin closed loop for 78 different initial conditions. 
In the second scenario we perturb the state by 1%, for 68 different initial conditions. We solve 
15 consecutive problems per initial condition and subsequently compare the average number 
of iterations as well as the generated hitting times per problem solve with and without warm¬ 
starting the dual variables. The results are summarized in FigureIt is evident that warm-strating 
consistently reduces the number of iterations in both cases. 

Finally, the active reduction in the horizon length thanks to warm-starting is illustrated in 
Figure for four different initial states in the case of the 1% perturbation. There is an (almost 
monotonic) decreasing trend, with the horizon finally shrinking to zero when the initial state 
resides in the positively invariant set. 


VII. Conclusion 

This work presented an algorithmic scheme capable of solving the constrained linear quadratic 
regulator problem in real time. The algorithm is an accelerated version of Forward-Backward 
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Box Plots: Cold VS Warm starting 



Fig. 6. ‘Box & Whiskers’ plot for the average number of iterations of the warm-starting policy, given 0.5% and 1% 
uniform perturbations of the initial state. The horizontal line inside the box corresponds to the median, the edges 
of the box are the 25*^ and 75*^ percentiles while the horizontal lines outside the boxes correspond to the most 
extreme data points not considered outliers. Finally, the colored dots correspond to mean values. One can observe 
that in the case of warm-starting the means are located outside the extreme points band, suggesting the existence 
of large outliers. However, warm-starting improves the performance of the suggested method in both cases, in all 
the depicted statistical measures. 


Hitting times evoiution for receding horizon control 



Fig. 7. Hitting times evolution for 50 subsequent solves, warm-started at the previous (shifted) dual optimizer. Four 
different initial conditions are depicted with the solid lines. The linear rate is depicted (for comparison) with the 
blue triangles. Note that the horizon decreases almost linearly, resulting from the relatively small perturbations of 
the initial state. The case depicted in red corresponds to an outlier where the perturbed initial condition rendered 
the previously computed horizon an infeasible option, leading to an increase of its length. 


Splitting, a popular proximal gradient scheme. The approach is to write the problem in its 
condensed form and dualize, which leads to the minimization of an infinite-dimensional quadratic 
functional subject to non-positivity constraints. The resulting infinite dimensional problem can 
be tackled in finite dimensions by observing that the dual sequence has always only finitely many 
non-zero elements. The proposed algorithm makes no use of terminal invariant sets and provably 
converges to the optimal solution of the infinite-horizon problem. Regarding the implementation 
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aspects, the algorithm can be highly competitive since it enjoys the convergence properties of 
optimal first order methods whose each iteration is computationally cheap. In addition, it requires 
minimal a priori information since the most crucial quantities are computed online, and there are 
no unreasonable or conservative assumptions on the problem’s structure. On the other hand, the 
number of iterations needed for convergence can be quite large for very unstable (large) systems. 
To this end, preconditioning of the problem has to be considered, with the obvious difficulty of 
computing preconditioners for infinite dimensional operators. Some steps towards this direction 
are taken in, e.g., p7| , but this problem remains a topic of future research. 

The algorithm used, AFBS, requires the objective function to be written as the sum of a 
smooth and a (possibly) nonsmooth function whose proximal operator is inexpensive to compute. 
Smoothness of the dual problem is recovered from strong convexity of the primal. The polytopic 
constraint set becomes a nonnegativity constraint when looking at the dual variables, resulting 
in a simple proximal operator. 

A direct extension of the existing scheme is the case of (some states) tracking a reference 
signal. By casting the tracking problem in a delta-formulation with respect to the inputs and states, 
the resulting problem retains the same structure as the regulation problem, with the difference 
of varying constraint sets due to the generated steady-state points. 

Appendix A 

Required Operator Theory 

For the sake of clarity we introduce some notation and definitions from operator theory. 
The subsequent results hold for general real Hilbert spaces, including the special case of 
we consider. For an in-depth treatment of (monotone) operator theory, the interested reader is 
referred to 0 and [ |28| . We write variables in normal font, and we use the bold font to describe 
the infinite-dimensional variables we are manipulating in our problem description. 

Definition 1: The /^-weighted (or real Hilbert space l-L is defined by 

H = |z = (^i)igpj ; ^ \\zi\\lw^ < ooj , w > 0 . 

Definition 2: A linear operator (mapping) F: Hi —> 7^2 between two Hilbert spaces is said 
to be bounded if the operator norm ||F|| of F, defined as 

||F|| := sup ||Fa:||^, , 

l|a;||«i=i 

satisfies ||F|| < oo. The set of bounded operators between two Hilbert spaces Hi and H 2 is 
denoted as -B(Hi,H 2 ). 

Theorem 2: Let Hi,H 2 be real Hilbert spaces and F G The adjoint of F is the 

unique operator F* G i3(H2,Hi) that satisfies 

{Fx,y) = {x,F^y) Vx G Hi, Vy G H 2 • 

Moreover, ||F|| = ||F*||. 

Subsequently, we introduce the notions of weak and strong convergence. 

Definition 3: Let H be a Hilbert space. We say that (x^) converges weakly to x if Vi/ G H 

( 1 /, x^) {y,x). We denote weak convergence as x^ ^ x. 
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Definition 4: Let be a sequenee in H. Then eonverges strongly to x if ||a;^ — 

x|| 0. We denote strong eonvergenee as x^ —>■ x. 

Definition 5: An operator F : 14, ^ 14 is positive-definite if it is bounded, F = F* and 
{Fx,x) > a||a:||-H for some a > 0, for all a: G "H . 

We reeall that a positive definite operator F is invertible and the inverse operator F~^ is bounded. 

Definition 6: The proximity operator of a proper lower semi-eontinuous eonvex funetion 
/; H ^ M is 

prox f{y) = argmin f{x) + -||a: - y\\^ 
x&n L ^ 



Appendix B 

Convergence oe the accelerated input-state sequences oe Algorithm!^ 

In this seetion we present the result that allows us to prove boundedness of the horizon 
sequenee T^, point (v) of Theorem in Seetion IV It is suffieient to show weak eonvergenee 
of the sequenee to the state optimizer x°°, whieh is required in the proof of point (v) of 
Theorem [TJ 

The theorem is stated at the end of the Appendix as a sequenee of intermediate results. We 
show that: 

1) The relaxed dual sequenee A eonverges weakly to a dual minimizer A°°. 

2) Provided that the operator C is bounded, the sequenee Hu{= C*X — h) eonverges 
weakly to Hu^, and eonverges weakly to u^. From strong duality, eonverges to the 
primal optimizer u°°. 

3) Weak eonvergenee of the aeeelerated state sequenee x^ to x°° follows direetly. 


Weak eonvergenee of the relaxed sequenee A follows from Corollary 2 of 112|, whieh states 
that the error sequenee (||A^ — eonverges to zero with rate l/k4. We state the result 

below. ^ 

Lemma 3: The relaxed sequenee A eonverges weakly to A°°. 

Proof: Sinee ||A^ — A^“^|p —0 and at is bounded we also have = \/c^(A^- 
Sinee strong eonvergenee implies weak eonvergenee we have that (iz^, y) 0, Vt/ G Hx. 


, k—l\ 


0 . 


The relaxed sequenee of duals A ean be written as A = A^ 
A^ ^ A“, we have that 


+ 




Consequently, sinee 


{X ,y) = {X’^ + Vc^iy^,y) 

= (A^, y) + y) X° 


for all y G Hx and henee A ^ A“. ■ 

Lemma 4: The sequenee eonverges weakly to u°°. 

Proof: Writing down the relation between ii and A from Lemma in terms of the operators, 
we get ii^ = H~^{C'^WX — h). Similarly we have = H~^{C*WX^ — h). Now, from 

n . 7- - - . - _ m . c k 


Theorem [ij we have —)• u°° and from Lemma we have A — A ^ 0. Therefore, sinee 
C*, W and H~^ are bounded operators (see Lemmas and in Appendix C) and sinee weak 
eonvergenee is preserved under bounded linear mappings, we eonelude that ii^ u°°. 


Lemma 5: The sequenee eonverges weakly to x 
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Proof: Exactly as we did at point (iv) of Theorem the aeeelerated state sequenee ean be 
written as 

= Axinit + . 

Weak eonvergenee u°° and boundedness of B prove weak eonvergenee of the aeeelerated 
state sequenee to x°°. ■ 


Appendix C 

Boundedness of several operators 

Lemma 6: The operators C and W are bounded. 

Proof: Boundedness of W is trivial sinee it is a diagonal operator with non-inereasing 
elements on the diagonal. 

The operator C ean be expressed as the following sum: 


C = 


■ a 

0 ■■■ ■ 


■ 0 

0 ■■■ ■ 

0 

0 ■■■ 


CxB 

0 

0 

Cu ■■■ 

+ 

0 

0 ■■■ 

0 

0 ■■■ 


CxAB 

CxB ■■■ 









and by using the triangle inequality, we have 

||C*|| < (Tmax(C'n) + \\Ca 

We thus have to show that y = C^u is bounded, i.e.. 


(19) 


sup \\y\\nx < oo- 

In order not to earry the zero rows of Cx, we define y G Tix by dropping the zero elements of y. 
This infinite-dimensional veetor eonsists of the elements ft = [Cx]iU = X]j=c) ^ 

Note that 

sup WvW'Hx = sup WyWn^ . 
l|aa||H„ = l l|aa||H„ = l 

Foeusing on the operator of interest, we have that 

\\y\\n^ = 


\ *=1 


w = 


i—1 


\ i=l j=0 


2 
Jll2 


( 20 ) 


Then 


oo i—1 

i=l j=0 

OO i—1 

i=l j=0 

OO i—1 

i=l j=o 


2 

7112 


.jW II2 


2 

7112 5 
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where we introdueed w = Uj = UjW^ with ii G and A = Aw. Observing the above 
expression, one ean identify that CxA^~^~^Buj is the convolution sum of the impulse 

response of the system 


S := 


A 

B 

a 

0 


with an input u. More speeifieally, borrowing the notation from [29| we denote the impulse 
response of S as = CxA’‘~^Bsi, where Sj = 1, i > 0 is the unit step funetion. Then 
the eonvolution operator is defined as the linear map S-^: ii —)■ y with |/j = * u)^ = 

X]!=o hs^i-jUj = (S'sfi)j. Thus we have that 


WWh^ = w 


oo 

2=1 


{S^u\\\l = w^S^u 


\\y\\nx =< w sup • 

II'“I|2<1 

From the definition of the indueed 2-norm of S, denoted here as ||S|| 2 , we have that sup ||5'sfi||2 = 
||S|| 2 ||'u ||2 and by assuming (without loss of generality) that ||fi ||2 = we end up 

having that <t?)||S|| 2 . Finally, the operator is bounded by the T-Loo norm of the transfer 

matrix H^{z), or < m)'Hoo(E). Subsequently, we have from (19) that the operator C is 

bounded by crmax(C'„) + wUoo{^)- ■ 

Remark 7: Following the diseussion from Seetion the weight w ean be ehosen to render 
any unstable system stable by shrinking the eigenvalues of the matrix A {A = wA). 

Lemma 7: The operator H~^ is bounded and ||iT|| < l/X m^A R). 

Proof: The operator H is given by iT = B*QB + R; see 0. Sinee B*QB is positive 
semidefinite (i.e., {B*QBu, u) > 0 for all u E 'Hu) and R is positive definite aeeording to 
Definition we eonelude that H is positive definite and henee has a bounded inverse. 

In order to eompute the bound for we observe that \\(B*QB + < 

1/Amin(-R), whieh eoneludes the proof. ■ 


Appendix D 

Backtracking stepsize ruee 


Theorem 1 in [ 12|, proves eonvergenee in the funetion values of F = f + g as well as weak 
eonvergenee of the iterates under the assumption of a fixed stepsize p = with ||V/(Ai) — 
V/(A2)|| < L(/)||Ai — A 2 II. As diseussed in Seetion |v[ an upper bound for the global Lipsehitz 
eonstant ean be overly eonservative and lead to a small stepsize. In this Appendix we briefly 
revise the baektraeking stepsize rule that allows for loeal estimates of the eurvature of / and show 


that the results of Theorem 1 also hold in this ease. The arguments are in line with [16|. Sinee 


the eonvergenee theory developed in [16| applies to general Hilbert spaees (see Remark 2.1), 
the following hold in the spaee of interest. 

We denote as /(■) the smooth and gf) the nonsmooth eonvex funetions of interest. Note 
that in our ease f = h* and g = For any L > 0 eonsider the quadratie approximation of 
F(A) = /(A) + g{X) at a point y: 


Ql{X, y) := f{y) + (A - z/, V/(z/)) + -||A - y\\'^ + ^(A) 


( 21 ) 
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We also define the unique minimizer parametrized by the point y as 
Piiv) := argmin{Qi(A,?/)} 


= argmin |^(A) + ^||A - - ^V/(?/)^ 


= proxg 1^?/ - -Vf{y) 

which is the basic step of Algorithmic i.e.. Step 3. The following holds: 
Lemma 8 (Lemma 2.3 Let y E LLx and L > 0 be such that 


F{pL{y)) < QL{PL{y),y) ■ 


Then for any A G T-Lx, 


F(A) - F(pL{y)) > -\\PL{y) - yf + L{y - X,PL{y) - y) ■ 


The backtracking procedure as described in [16| is as follows: 


( 22 ) 

(23) 

(24) 


Algorithm 3 Backtracking for stepsize computation 

0: Take > 0, some p > I, and A° G Hx- 

repeat 

1: Find the smallest nonnegative integer such that with L = p" L^~^ 
until Fipiiy’^)) < Qi(pi(y),y) 


Note that for any L^ = L generated by Algorithmic Lemma holds. 

We also have the following instrumental Lemma from p^ : 

Lemma 9 (Lemma 1 4^2| /).' Let L > L(f), X,y eHx and piiy) := pro^g{y - ^Vf{y)). 
Then for all A 

F{PL{y)) + ^\\pL{y) - yV < F{y) + - i/f • 

Lemma is crucial for proving point (i) of Theorem |C More specifically, supposing that 
Lemma|Cholds, convergence of the function values with rate l//c^ can be proven under no further 
assumptions using Theorem 2 and Corollary 1 in p^ . We will not repeat the aforementioned 
theorems here due to space limitations. What we are going to show instead is that all stepsizes 
= 1/L^ generated from Algorithm |C satisfy Lemma |C 

Lemma 10: Consider F = f+g, with f = h* and p = (5_ as defined in Section |nj The iterates 
A^ generated from Algorithm |C with a backtracking stepsize rule generated from Algorithm [C 
satisfy: 

Proof: All local Lipschitz estimates L generated from Algorithm |C satisfy F(pi(y^)) < 
QL{PL{y)iy) and, consequently, from Lemma |C 


F(A) - F(pi(y)) > ^Wpiiy) - y\? + L{y - X,pi(y) - y) , 
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2 

J- (F(A) - F{pi{y))) > Wpiiy) - yf + 2{y - X^piiy) - y) ■ 

It follows from the Pythagorean theorem that 

||6 — a|p + 2(6 — a, a — c) = ||6 — c|l^ — ||a — c||^ , 

and henee 

j- (F(A) - F{pi{y))) > \\pi{y) - yf - \\y - Af , 

from whieh Lemma [^follows. Theorem 2 and Corollary 1 of [ 12| then lead to the desired result. 
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